#include <mex.h>
#include <math.h>

void ED_Bandfun1(double *dCdt, double *dC1dt, double *dC2dt, double t, double *ffd,
                double *C, double *C1, double *C2, int NE, int NP,
                double cf, double cf1, double cf2, double *dpp, double *AZ, double *BZ)
{
double L, d, dp, et, ep, eb, F, VS, u, W, DL, DL1, DL2, RI, RI1, RI2; 
int    i, ii, k;
double Q, S, K1, K2, z1, z2;

double *lambdaS, *lambda1S, *lambda11, *lambda2S, *lambda22;

double *D, *D1, *D2, *CC, *CC1, *CC2, *DD, *DD1, *DD2, *x, *x1, *x2;
double s = 0.0, ss = 0.0, s1 = 0.0, s2 = 0.0, s3 = 0.0, s4 = 0.0;
double Le;
int    M, lgt1;

lambdaS  = mxCalloc(NE*NP, sizeof(double));
lambda1S = mxCalloc(NE*NP, sizeof(double));
lambda11 = mxCalloc(NE*NP, sizeof(double));
lambda2S = mxCalloc(NE*NP, sizeof(double));
lambda22 = mxCalloc(NE*NP, sizeof(double));
D  = mxCalloc(NE+1, sizeof(double));
D1 = mxCalloc(NE+1, sizeof(double));
D2 = mxCalloc(NE+1, sizeof(double));
CC  = mxCalloc(NE+1, sizeof(double));
CC1 = mxCalloc(NE+1, sizeof(double));
CC2 = mxCalloc(NE+1, sizeof(double));
DD  = mxCalloc(NE+1, sizeof(double));
DD1 = mxCalloc(NE+1, sizeof(double));
DD2 = mxCalloc(NE+1, sizeof(double));
x  = mxCalloc(NE+1, sizeof(double));
x1 = mxCalloc(NE+1, sizeof(double));
x2 = mxCalloc(NE+1, sizeof(double));

et = dpp[0]; ep = dpp[1]; L = dpp[2]; d = dpp[3]; dp = dpp[4]; W = dpp[5]; //VS = dpp[6];
Q = ffd[0]; S = ffd[1]; K1 = ffd[2]; K2 = ffd[3]; z1 = ffd[4]; z2 = ffd[5]; VS = ffd[6]; DL1 = ffd[7]; DL2 = ffd[8];

eb = (et-ep)/(1-ep);
u = 4.0*VS/(d*d*3.14159*eb);/*Stromgeschwindigkeit*/
//F = (1-et)/et;
F = (1-eb)/eb;
DL = (dp*u*eb)/(0.2+0.011*pow(dp*u/W,0.48));/*alxiale Diffusionskoeffizient*/
Le = 1.0/NE;
lgt1 = NE*NP-1;
M = NP+2;
//printf("DL = %f, DL1 = %f, DL2 = %f\n", DL, DL1, DL2);

for (i=0; i<=lgt1; i++)
{
    dCdt [i] = 0.0; /*Salz*/
    dC1dt[i] = 0.0; /*Protein 1*/
    dC2dt[i] = 0.0; /*Protein 2*/
}

for (i=1; i<=NP; i++)
{
    s  = s +AZ[M*i]*C [i-1]; /*Salz*/
    s1 = s1+AZ[M*i]*C1[i-1]; /*Protein 1*/
    s2 = s2+AZ[M*i]*C2[i-1]; /*Protein 2*/
    ss = ss+AZ[NP+1+M*i]*C [(NE-1)*NP+i-1]; /*Salz*/
    s3 = s3+AZ[NP+1+M*i]*C1[(NE-1)*NP+i-1]; /*Protein 1*/
    s4 = s4+AZ[NP+1+M*i]*C2[(NE-1)*NP+i-1]; /*Protein 2*/
}
D [0] = -u*L*Le*cf -DL*s ; /*Salz*/
D1[0] = -u*L*Le*cf1-DL1*s1; /*Protein 1*/
D2[0] = -u*L*Le*cf2-DL2*s2; /*Protein 2*/
D [NE] = -ss; /*Salz*/
D1[NE] = -s3; /*Protein 1*/
D2[NE] = -s4; /*Protein 2*/
CC [0] = DL *AZ[M*(NP+1)]/(DL *AZ[0]-u*L*Le);
CC1[0] = DL1*AZ[M*(NP+1)]/(DL1*AZ[0]-u*L*Le);
CC2[0] = DL2*AZ[M*(NP+1)]/(DL2*AZ[0]-u*L*Le);
DD [0] = D [0]/(DL *AZ[0]-u*L*Le); /*Salz*/
DD1[0] = D1[0]/(DL1*AZ[0]-u*L*Le); /*Protein 1*/
DD2[0] = D2[0]/(DL2*AZ[0]-u*L*Le); /*Protein 2*/
// D [0] = -u*L*Le*cf -eb*DL*s ; /*Salz*/
// D1[0] = -u*L*Le*cf1-eb*DL*s1; /*Protein 1*/
// D2[0] = -u*L*Le*cf2-eb*DL*s2; /*Protein 2*/
// D [NE] = -ss; /*Salz*/
// D1[NE] = -s3; /*Protein 1*/
// D2[NE] = -s4; /*Protein 2*/
// CC [0] = eb*DL*AZ[M*(NP+1)]/(eb*DL*AZ[0]-u*L*Le);
// CC1[0] = eb*DL*AZ[M*(NP+1)]/(eb*DL*AZ[0]-u*L*Le);
// CC2[0] = eb*DL*AZ[M*(NP+1)]/(eb*DL*AZ[0]-u*L*Le);
// DD [0] = D [0]/(eb*DL*AZ[0]-u*L*Le); /*Salz*/
// DD1[0] = D1[0]/(eb*DL*AZ[0]-u*L*Le); /*Protein 1*/
// DD2[0] = D2[0]/(eb*DL*AZ[0]-u*L*Le); /*Protein 2*/

for(k=1; k<=NE-1; k++)
{
    s = 0.0; ss = 0.0; s1 = 0.0; s2 = 0.0; s3 = 0.0; s4 = 0.0;
    for(i=1; i<=NP; i++)
    {
        s  = s +AZ[M*i]*C [k*NP+i-1];   /*Salz*/
        ss = ss+AZ[NP+1+M*i]*C [(k-1)*NP+i-1]; /*Salz*/
        s1 = s1+AZ[M*i]*C1[k*NP+i-1]; /*Protein 1*/
        s2 = s2+AZ[NP+1+M*i]*C1[(k-1)*NP+i-1]; /*Protein 1*/
        s3 = s3+AZ[M*i]*C2[k*NP+i-1];
        s4 = s4+AZ[NP+1+M*i]*C2[(k-1)*NP+i-1];
    }
    D [k] = s -ss; /*Salz*/
    D1[k] = s1-s2; /*Protein 1*/
    D2[k] = s3-s4; /*Protein 2*/
    CC [k] = -AZ[M*(NP+1)]/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC [k-1]);
    CC1[k] = -AZ[M*(NP+1)]/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC1[k-1]);
    CC2[k] = -AZ[M*(NP+1)]/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC2[k-1]);
    DD [k] = (D [k]-AZ[NP+1]*DD [k-1])/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC [k-1]); /*Salz*/
    DD1[k] = (D1[k]-AZ[NP+1]*DD1[k-1])/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC1[k-1]); /*Protein 1*/
    DD2[k] = (D2[k]-AZ[NP+1]*DD2[k-1])/(AZ[NP+1+M*(NP+1)]-AZ[0]-AZ[NP+1+M*0]*CC2[k-1]); /*Protein 2*/
}
DD [NE] = (D [NE]-AZ[NP+1]*DD [NE-1])/(AZ[NP+1+M*(NP+1)]-AZ[NP+1]*CC [NE-1]); /*Salz*/
DD1[NE] = (D1[NE]-AZ[NP+1]*DD1[NE-1])/(AZ[NP+1+M*(NP+1)]-AZ[NP+1]*CC1[NE-1]); /*Protein 1*/
DD2[NE] = (D2[NE]-AZ[NP+1]*DD2[NE-1])/(AZ[NP+1+M*(NP+1)]-AZ[NP+1]*CC2[NE-1]); /*Protein 2*/
x [NE] = DD [NE]; /*Salz*/
x1[NE] = DD1[NE]; /*Protein 1*/
x2[NE] = DD2[NE]; /*Protein 2*/

for(k=NE; k>=1; k--)
{
    x [k-1] = DD [k-1]-CC [k-1]*x [k]; /*Salz*/
    x1[k-1] = DD1[k-1]-CC1[k-1]*x1[k]; /*Protein 1*/
    x2[k-1] = DD2[k-1]-CC2[k-1]*x2[k]; /*Protein 2*/
    /*printf("x [%d] = %f\n",k, x [k]);
    printf("x1[%d] = %f\n",k, x1[k]);
    printf("x2[%d] = %f\n",k, x2[k]);*/
}

for(i=0; i<=lgt1; i++) 
{
   lambdaS [i] = 2*S*C[i]*pow(4*S*C[i]*C[i],-0.5); /*Salz*/
   lambda1S[i] = -z1*K1*pow(C[i],-z1-1); /*Protein 1*/
   lambda11[i] = K1*pow(C[i],-z1); /*Protein 1*/
   lambda2S[i] = -z2*K2*pow(C[i],-z2-1); /*Protein 2*/
   lambda22[i] = K2*pow(C[i],-z2); /*Protein 2*

/*printf("C [%d] = %f\n",i,C [i]);
printf("C1[%d] = %f\n",i,C1[i]);
printf("C2[%d] = %f\n",i,C2[i]);
printf("lambdaS [%d] = %f\n",i,lambdaS [i]);
printf("lambda1S[%d] = %f\n",i,lambda1S[i]);
printf("lambda11[%d] = %f\n",i,lambda11[i]);
printf("lambda2S[%d] = %f\n",i,lambda2S[i]);
printf("lambda22[%d] = %f\n",i,lambda22[i]);*/
}
for(k=1; k<=NE; k++)
{
   for(i=2; i<=NP+1; i++)
   {   
      s = 0.0; ss = 0.0; s1 = 0.0; s2 = 0.0; s3 = 0.0; s4 = 0.0;
      for(ii=1; ii<=NP; ii++)
      {
         s  = s +BZ[i-1+M*ii]*C [(k-1)*NP+ii-1]; /*Salz*/
         ss = ss+AZ[i-1+M*ii]*C [(k-1)*NP+ii-1]; /*Salz*/
         s1 = s1+BZ[i-1+M*ii]*C1[(k-1)*NP+ii-1]; /*Protein 1*/ 
         s2 = s2+AZ[i-1+M*ii]*C1[(k-1)*NP+ii-1]; /*Protein 1*/
         s3 = s3+BZ[i-1+M*ii]*C2[(k-1)*NP+ii-1]; /*Protein 2*/
         s4 = s4+AZ[i-1+M*ii]*C2[(k-1)*NP+ii-1]; /*Protein 2*/
      }
      RI  = DL *(BZ[i-1]*x [k-1]+s +BZ[i-1+M*(NP+1)]*x [k])/(L*L*Le*Le)-u*(AZ[i-1]*x [k-1]+ss+AZ[i-1+M*(NP+1)]*x [k])/(L*Le); /*Salz*/
      RI1 = DL1*(BZ[i-1]*x1[k-1]+s1+BZ[i-1+M*(NP+1)]*x1[k])/(L*L*Le*Le)-u*(AZ[i-1]*x1[k-1]+s2+AZ[i-1+M*(NP+1)]*x1[k])/(L*Le); /*Protein 1*/
      RI2 = DL2*(BZ[i-1]*x2[k-1]+s3+BZ[i-1+M*(NP+1)]*x2[k])/(L*L*Le*Le)-u*(AZ[i-1]*x2[k-1]+s4+AZ[i-1+M*(NP+1)]*x2[k])/(L*Le); /*Protein 2*/
      dCdt [(k-1)*NP+i-2] = RI/(1+F*lambdaS[(k-1)*NP+i-2]); /*Salz*/
      dC1dt[(k-1)*NP+i-2] = (RI1-F*lambda1S[(k-1)*NP+i-2]*RI/(1+F*lambdaS[(k-1)*NP+i-2]))/(1+F*lambda11[(k-1)*NP+i-2]); /*Protein 1*/
      dC2dt[(k-1)*NP+i-2] = (RI2-F*lambda2S[(k-1)*NP+i-2]*RI/(1+F*lambdaS[(k-1)*NP+i-2]))/(1+F*lambda22[(k-1)*NP+i-2]); /*Protein 2*/
//       RI  = eb*DL*(BZ[i-1]*x [k-1]+s +BZ[i-1+M*(NP+1)]*x [k])/(L*L*Le*Le*et)-eb*u*(AZ[i-1]*x [k-1]+ss+AZ[i-1+M*(NP+1)]*x [k])/(L*Le*et); /*Salz*/
//       RI1 = eb*DL*(BZ[i-1]*x1[k-1]+s1+BZ[i-1+M*(NP+1)]*x1[k])/(L*L*Le*Le*et)-eb*u*(AZ[i-1]*x1[k-1]+s2+AZ[i-1+M*(NP+1)]*x1[k])/(L*Le*et); /*Protein 1*/
//       RI2 = eb*DL*(BZ[i-1]*x2[k-1]+s3+BZ[i-1+M*(NP+1)]*x2[k])/(L*L*Le*Le*et)-eb*u*(AZ[i-1]*x2[k-1]+s4+AZ[i-1+M*(NP+1)]*x2[k])/(L*Le*et); /*Protein 2*/
//       dCdt [(k-1)*NP+i-2] = RI/(1+F*lambdaS[(k-1)*NP+i-2]); /*Salz*/
//       dC1dt[(k-1)*NP+i-2] = (RI1-F*lambda1S[(k-1)*NP+i-2]*RI/(1+F*lambdaS[(k-1)*NP+i-2]))/(1+F*lambda11[(k-1)*NP+i-2]); /*Protein 1*/
//       dC2dt[(k-1)*NP+i-2] = (RI2-F*lambda2S[(k-1)*NP+i-2]*RI/(1+F*lambdaS[(k-1)*NP+i-2]))/(1+F*lambda22[(k-1)*NP+i-2]); /*Protein 2*/
      /*printf("dCdt[(%d-1)*NP+%d-2] = %f\n",k,i, dCdt[(k-1)*NP+i-2]);
      printf("dC1dt[(%d-1)*NP+%d-2] = %f\n",k,i, dC1dt[(k-1)*NP+i-2]);
      printf("dC2dt[(%d-1)*NP+%d-2] = %f\n",k,i, dC2dt[(k-1)*NP+i-2]);*/
    }
}  
}
    
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double t, cf, cf1, cf2;
    double *dpp, *ffd;
    double *dCdt, *dC1dt, *dC2dt;
    double *C, *C1, *C2;
    double *AZ, *BZ;
    int NE, NP;
    t = mxGetScalar(prhs[0]);
    ffd = mxGetPr(prhs[1]);
    C  = mxGetPr(prhs[2]);
    C1 = mxGetPr(prhs[3]);
    C2 = mxGetPr(prhs[4]);
    NE = mxGetScalar(prhs[5]);
    NP = mxGetScalar(prhs[6]);
    cf  = mxGetScalar(prhs[7]);
    cf1 = mxGetScalar(prhs[8]);
    cf2 = mxGetScalar(prhs[9]);
    dpp = mxGetPr(prhs[10]);
    AZ = mxGetPr(prhs[11]);
    BZ = mxGetPr(prhs[12]);
    plhs[0] = mxCreateDoubleMatrix(NE*NP, 1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(NE*NP, 1, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(NE*NP, 1, mxREAL);
    dCdt  = mxGetPr(plhs[0]);
    dC1dt = mxGetPr(plhs[1]);
    dC2dt = mxGetPr(plhs[2]);
    ED_Bandfun1(dCdt, dC1dt, dC2dt, t, ffd, C, C1, C2, NE, NP, cf, cf1, cf2, dpp, AZ, BZ);
}